General Plots for M.1. (Autobiographical Memories)

Initial Dataset

1 Read

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Read
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

## Read files
file <- "Date Complete M1 v.13 siPPGGSRamilaza.sav"

# setwd(folder)
Data <- rio::import(file)

2 Make data frame

Data %>%
  dplyr::select(-Nume) %>%
    DT::datatable(
      extensions = 'Buttons',
      options = list(pageLength = 10,
                     scrollX='500px',
                     dom = 'Bfrtip',
                     buttons = c('excel', "csv")))

2.1 Exclude Protocol 8 (mother)

Data <- 
  Data %>%
  filter(P != 8)

3 Define Functions

## Func t test si boxplot simplu
func_t_box <- function(df, ind, pre_var, post_var, facet = FALSE, xlab = ""){  
  if(facet){
    facet <- "Protocol"
  }else{
    facet <- NULL
  }
  
  df_modif <-
    df %>%
    select(ind, P, pre_var, post_var) %>% 
    tidyr::drop_na() %>%
    gather(pre_var, post_var, key = "PrePost", value = "value") %>% 
    mutate_at(vars(c(1, 2)), funs(as.factor)) %>% 
    mutate(PrePost = factor(PrePost, levels = c(pre_var, post_var))) 
  
  if(!is.null(facet)){
    df_modif <-
      df_modif %>%
      group_by(P) %>%
      mutate(Protocol = paste0("Protocol = ", P, ", n = ", n()))
  }
  
  stat_comp <-
    df_modif %>% 
    do(tidy(t.test(.$value ~ .$PrePost,
                   paired = TRUE,
                   data=.)))
  
  plot <- 
    ggpubr::ggpaired(df_modif, x = "PrePost", y = "value", id = ind, 
                     color = "PrePost", line.color = "gray", line.size = 0.4,
                     palette = c("#00AFBB", "#FC4E07"), legend = "none",
                     facet.by = facet, ncol = 3, 
                     xlab = xlab) +
    stat_summary(fun.data = mean_se,  colour = "darkred") +
    ggpubr::stat_compare_means(method = "t.test", paired = TRUE, label.x = as.numeric(df_modif$PrePost)-0.4, label.y = max(df_modif$value)+1) + 
    ggpubr::stat_compare_means(method = "t.test", paired = TRUE, label = "p.signif", comparisons = list(c(pre_var, post_var)))
  
  print(stat_comp)
  cat("\n")                      
  print(plot)
  cat("\n")
  plot.new()                     # Need this workaround for interleaving tables and plots in R Markdown, within loop
  dev.off()
}

heat_cor_plotly <- function(df, x_vars = NULL, y_vars = NULL, low_color = "cyan",  high_color = "red",  ...){   
  # inherit type = c("pearson","spearman") from Hmisc::rcorr() 
  library(ggplot2)
  library(plotly)
  library(reshape2)
  library(Hmisc)
  
  # use all numeric columns only, print message if non-numeric are found
  numeric_cols <- unlist(lapply(df, is.numeric))
  if(!all(numeric_cols)) message("Warning: Non-numeric columns were excluded!")
  df <- df[, numeric_cols]
  
  df_mat <- as.matrix(df)
  rt <- Hmisc::rcorr(df_mat, ...)
  
  # extract correlations, p-values and merge into another dataframe
  mtlr <- reshape2::melt(rt$r, value.name = "Correlation")
  mtlp <- reshape2::melt(rt$P, value.name = "P-Value")
  
  mtl <- merge(mtlr, mtlp)
  
  # give possibility to prune the correlation matrix
  if(!is.null(x_vars)){
    mtl <- mtl[(mtl$Var1 %in% x_vars), ]
  }
  if(!is.null(x_vars)){
    mtl <- mtl[(mtl$Var2 %in% y_vars), ]
  }
  
  # want to avoid scientific notetion, but this doesnt work as numeric
  # mtl$Correlation <- as.numeric(format(mtl$Correlation, digits = 4, scientific = FALSE))  # doesnt work
  # mtl$`P-Value` <- as.numeric(format(mtl$`P-Value`, digits = 4, scientific = FALSE)) 
  options(scipen = 999)
  mtl$Correlation <- round(mtl$Correlation, 3)
  mtl$`P-Value` <- round(mtl$`P-Value`, 3)

  gx <-
    ggplot2::ggplot(mtl, 
           aes(Var1, Var2, 
               fill = Correlation,  
               text = paste("P-val = ", `P-Value`))) +
    ggplot2::geom_tile() + 
    ggplot2::scale_fill_gradient(low = low_color,  high = high_color, limits = c(-1, 1), breaks = c(-1, -.5, 0, .5, 1)) +
    ggplot2::theme_minimal() +
    {if(any(nchar(names(df)) > 6)) ggplot2::theme(axis.text.x = element_text(angle = 90, hjust = 1))}  # vertical x axis labels if lenghty
  plotly::ggplotly(gx)  
}

4 Plot Age

## Dodged Bar plot of Age and Gender
Data  %>%
  mutate(Varta_categ = cut(Varsta, 
                           breaks=c(-Inf, 25, 30, 35, 40, 45, 50, 55, 60, Inf), 
                           labels=c("<25","25-29","30-34", "35-39", "40-44", "45-49", "50-54", "55-59", "60>"), 
                           right = FALSE)) %>%  
  mutate(Varsta = as.factor(Varsta),
         Gen = as.factor(as.character(Gen))) %>%
  mutate(Gen = forcats::fct_recode(Gen, "femin" = "1", "masculin" = "2")) %>%
  dplyr::count(Varta_categ, Gen, .drop = FALSE) %>%         # Group by, then count number in each group (dont drop 0 counts)
  mutate(pct = prop.table(n)) %>%                           # Calculate percent within each var
    ggplot(aes(x = Varta_categ, y = pct, fill = Gen, label = scales::percent(pct))) + 
      geom_col(position = position_dodge(preserve = "single"), stat = "identity",) +    # Don't drop zero count
      geom_text(position = position_dodge(width = .9),      # move to center of bars
                vjust = -0.5,                               # nudge above top of bar
                size = 3) + 
      scale_y_continuous(labels = scales::percent) +
      ggtitle("") +
      xlab("Varsta") + ylab("Percentage %") + 
      guides(fill = guide_legend(title = "Gen", ncol = 1)) + 
      scale_fill_grey(start = 0.8, end = 0.2, na.value = "red", aesthetics = "fill") +
      theme(legend.position = "right", legend.direction = "vertical", 
            legend.justification = c(0, 1), panel.border = element_rect(fill = NA, colour = "black"))

4.1 By Protocol

## Dodged Bar plot of Age and Gender by Protocol
Data  %>%
  mutate(Varta_categ = cut(Varsta, 
                           breaks=c(-Inf, 25, 30, 35, 40, 45, 50, 55, 60, Inf), 
                           labels=c("<25","25-29","30-34", "35-39", "40-44", "45-49", "50-54", "55-59", "60>"), 
                           right = FALSE)) %>%  
  mutate(Varsta = as.factor(Varsta),
         Gen = as.factor(as.character(Gen))) %>%
  mutate(Gen = forcats::fct_recode(Gen, "femin" = "1", "masculin" = "2")) %>%
  group_by(P) %>%
  dplyr::count(Varta_categ, Gen, .drop = FALSE) %>%         # Group by, then count number in each group (dont drop 0 counts)
  mutate(pct = prop.table(n)) %>%                           # Calculate percent within each var
    ggplot(aes(x = Varta_categ, y = pct, fill = Gen, label = scales::percent(pct))) +
      facet_wrap(~P, scales = "free", ncol = 1) +
      geom_col(position = position_dodge(preserve = "single"), stat = "identity",) +    # Don't drop zero count
      geom_text(position = position_dodge(width = .9),      # move to center of bars
                vjust = -0.5,                               # nudge above top of bar
                size = 3) + 
      scale_y_continuous(labels = scales::percent) +
      ggtitle("") +
      xlab("Varsta") + ylab("Percentage %") + 
      guides(fill = guide_legend(title = "Gen", ncol = 1)) + 
      scale_fill_grey(start = 0.8, end = 0.2, na.value = "red", aesthetics = "fill") +
      theme(legend.position = "right", legend.direction = "vertical", 
            legend.justification = c(0, 1), panel.border = element_rect(fill = NA, colour = "black"))

## Pie chart
Data  %>%
  mutate(Gen = as.factor(as.character(Gen))) %>%
  mutate(Gen = forcats::fct_recode(Gen, "femin" = "1", "masculin" = "2")) %>%
  group_by(Gen) %>%
  dplyr::summarise(counts = n()) %>%
  mutate(prop = round(counts*100/sum(counts), 1),
         lab.ypos = cumsum(prop) - .5*prop,
         Percent = paste0(prop, " %")) %>% 
  ggpubr::ggpie(x = "prop", label = "Percent",
                fill = "Gen", color = "white", 
                lab.pos = "in", lab.font = list(color = "white"),
                palette = "grey")

5 Analyses

5.1 Simple before-after analyses with t test

## Simple before-after analyses with t test
cat("#### VAS Stress")

5.1.0.1 VAS Stress

func_t_box(Data, ind = "ID", "Stres_pre", "Stres_post", facet = FALSE) 

NANA

null device 1

## Simple before-after analyses with t test
cat("#### VAS Stress")

5.1.0.2 VAS Stress

func_t_box(Data, ind = "ID", "Stres_pre", "Stres_post", facet = TRUE) 

NANA

null device 1

5.2 Correlations: Anotimpuri - Personality (without P6, P7)

dateplot_anopers <- Data[, c("P", "Primavara", "Vara", "Toamna", "Iarna")]
dateplot_anopers <- cbind(dateplot_anopers, Data[, 87:121])
dateplot_anopers <- subset(dateplot_anopers, P!=6 & P!=7)

COR <- Hmisc::rcorr(as.matrix(dateplot_anopers[, -1]))   
M <- COR$r[1:4, ]
P_MAT <- COR$P[1:4, ]
corrplot::corrplot(M, type = "upper", p.mat = P_MAT, sig.level = 0.05, insig = "blank", tl.col = "black", tl.cex = .7, cl.pos = "b", tl.srt = 45)

Error in data.frame(…, check.names = FALSE) : arguments imply differing number of rows: 150, 146

5.3 Correlations: Anotimpuri - Calitate Amintiri (without P6, P7)

dateplot1 <- Data[, c("P", "Primavara", "Vara", "Toamna", "Iarna", "Media_s1", "Media_s2", "Media_s3",  "SocDih_Part",  "SocDih_FamN",  "SocDih_FamInd",  "SocDih_Priet",  "SocDih_Amici",  "SocDih_Necun",  "SocDih_Antag",  "SocDih_TotAprop",  "SocDih_TotNeaprop", "STAI_T")] 
names(dateplot1) <- c("P", "Primavara", "Vara", "Toamna", "Iarna", "S1- Valenta", "S2 - Vividness", "S3 - Relevanta",  "Partener",  "Familie nucleu",  "Familie extinsa",  "Prieteni",  "Amici",  "Necunoscuti",  "Antagonisti",  "Toti Apropiatii",  "Toti Neapropiatii", "STAI_T")
dateplot1 <- subset(dateplot1, P!=6 & P!=7)

COR <- Hmisc::rcorr(as.matrix(dateplot1[,-1]))   
M <- COR$r
P_MAT <- COR$P
corrplot::corrplot(M, method = "number", type = "upper", p.mat = P_MAT, sig.level = 0.05, insig = "blank", tl.col = "black", tl.cex = .9, tl.srt = 45)  

Error in data.frame(…, check.names = FALSE) : arguments imply differing number of rows: 153, 136

heat_cor_plotly(dateplot1[,-1])

5.4 Correlations: Personality - Qualities of Memories (without P6, P7)

dateplot2 <- Data[, c(24, 40, 56, 87:121, 126)] 
names(dateplot2)[1:3] <- c("S1- Valenta", "S2 - Vividness", "S3 - Relevanta")

COR <- Hmisc::rcorr(as.matrix(dateplot2))   
M <- COR$r
P_MAT <- COR$P
corrplot::corrplot(M, type = "upper", p.mat = P_MAT, sig.level = 0.05, insig = "blank", tl.col = "black", tl.cex = .7, cl.pos = "b", tl.srt = 45)

Error in data.frame(…, check.names = FALSE) : arguments imply differing number of rows: 780, 741

heat_cor_plotly(dateplot2, x_vars = names(dateplot2)[1:3], y_vars = names(dateplot2)[-(1:3)])

5.5 Correlations: Social - Personality

dateplot3 <- Data[, c(131:139, 87:121)]
names(dateplot3)[1:9] <- c("Partener",  "Familie nucleu",  "Familie extinsa",  "Prieteni",  "Amici",  "Necunoscuti",  "Antagonisti",  "Toti Apropiatii",  "Toti Neapropiatii")

COR <- Hmisc::rcorr(as.matrix(dateplot3))   
M <- COR$r
P_MAT <- COR$P
corrplot::corrplot(M, type = "upper", p.mat = P_MAT, sig.level = 0.05, insig = "blank", tl.col = "black", tl.cex = .7, cl.pos = "b", tl.srt = 45)

Error in data.frame(…, check.names = FALSE) : arguments imply differing number of rows: 990, 946

heat_cor_plotly(dateplot3, x_vars = names(dateplot3)[1:9], y_vars = names(dateplot3)[-(1:9)])

6 Social

# names(Data[str_detect(colnames(Data), fixed("SocDih", ignore_case=TRUE))])
Data_soc <-
  Data %>%
  dplyr::select("ID", "P", starts_with("SocDih")) %>%
  haven::zap_formats(.) %>%                # to not get warning  
  haven::zap_labels(.) %>%                 # "attributes are not identical across measure variables"
  haven::zap_widths(.) %>%                 # on gather() because of SPSS
  dplyr::rename_all(list(~stringr::str_replace(., "SocDih_", ""))) %>% 
  gather(key = "Variable", value = "Value", -ID, -P) 
  
# Create a custom color scale for all ASCQ graphs
library(RColorBrewer)
myColors <- brewer.pal(10,"Set1")
names(myColors) <- levels(Data_soc$Variable)
colScale <- scale_colour_manual(name = "Variable", values = myColors)

# Plot
ggpubr::ggviolin(data = Data_soc, x = "Variable", y = "Value", fill = "Variable",
      add = "boxplot", add.params = list(fill = "white"),
      xlab = "", legend = "none") +
  colScale +                                                   # color scale here keep consistency of color with factor level
  stat_summary(fun.data = mean_se,  colour = "darkred") +         
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

# Plot faceted by Protocol
p <- 
  ggpubr::ggviolin(data = Data_soc, x = "Variable", y = "Value", fill = "Variable",
        add = "boxplot", add.params = list(fill = "white"),
        xlab = "", legend = "none") +
         
    colScale +                                                   # color scale here keep consistency of color with factor level
    stat_summary(fun.data = mean_se,  colour = "darkred") +         
    theme(axis.text.x = element_text(angle = 90, hjust = 1))

ggpubr::facet(p, facet.by = "P", ncol = 1, scales = "free_x")

6.1 Social - this doesnt make sense

# PerformanceAnalytics::chart.Correlation(Data[, c(8, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140)])

# Example of Simpsons Paradox
coplot(DifStres ~ SocDih_Amici | Media_s1,
       data = Data,
       rows = 1,
     panel = function(x, y, ...) {
          panel.smooth(x, y, span = .8, iter = 5,...)
          abline(lm(y ~ x), col = "blue")})

7 Varsta Amint - P1,P2,P3

Data_P1P2P3 <- 
  Data %>%
  filter(P %in% c("1", "2", "3")) %>%
  mutate(Med_amintvarsta = as.numeric(as.character(Med_amintvarsta)),
         Dif_Med_amintvarsta = Varsta - Med_amintvarsta)

PerformanceAnalytics::chart.Correlation(Data_P1P2P3[, c(8, 85, 122, 148)])

coplot(DifStres ~ Dif_Med_amintvarsta | Media_s1,
       data = Data_P1P2P3,
       rows = 1,
     panel = function(x, y, ...) {
          panel.smooth(x, y, span = .8, iter = 5,...)
          abline(lm(y ~ x), col = "blue")})

8 Protocol 3 - Social

Data_P3 <- 
  Data %>%
  filter(P == "3")

PerformanceAnalytics::chart.Correlation(Data_P3[, c(8, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140)])

# coplot(SocDih_Amici ~ DifStres | Media_s1, 
#        data = Data,
#        columns = 3,
#      panel = function(x, y, ...) {
#           panel.smooth(x, y, span = .8, iter = 5,...)
#           abline(lm(y ~ x), col = "blue") } )

9 Protocol 3 - Varsta Amint

Data_P3 <- 
  Data %>%
  filter(P == "3") %>%
  mutate(Med_amintvarsta = as.numeric(as.character(Med_amintvarsta)),
         Dif_Med_amintvarsta = Varsta - Med_amintvarsta)

PerformanceAnalytics::chart.Correlation(Data_P3[, c(8, 85, 122, 148)])

10 Session Info

R version 4.2.2 (2022-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 8.1 x64 (build 9600)

Matrix products: default

locale:
[1] LC_COLLATE=Romanian_Romania.1252  LC_CTYPE=Romanian_Romania.1252    LC_MONETARY=Romanian_Romania.1252
[4] LC_NUMERIC=C                      LC_TIME=Romanian_Romania.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] RColorBrewer_1.1-3         Hmisc_5.0-1                reshape2_1.4.4             plotly_4.10.1             
 [5] rio_0.5.29                 scales_1.2.1               ggpubr_0.6.0               rstatix_0.7.2             
 [9] broom_1.0.4                PerformanceAnalytics_2.0.4 xts_0.13.0                 zoo_1.8-11                
[13] psych_2.3.3                plyr_1.8.8                 lubridate_1.9.2            forcats_1.0.0             
[17] stringr_1.5.0              dplyr_1.1.1                purrr_1.0.1                readr_2.1.4               
[21] tidyr_1.3.0                tibble_3.2.1               ggplot2_3.4.2              tidyverse_2.0.0           
[25] papaja_0.1.1               tinylabels_0.2.3           pacman_0.5.3              

loaded via a namespace (and not attached):
 [1] colorspace_2.1-0    ggsignif_0.6.4      ellipsis_0.3.2      rprojroot_2.0.3     estimability_1.4.1  htmlTable_2.4.1    
 [7] parameters_0.20.3   base64enc_0.1-3     fs_1.6.1            rstudioapi_0.14     farver_2.1.1        DT_0.27            
[13] fansi_1.0.4         mvtnorm_1.1-3       mnormt_2.1.1        cachem_1.0.6        knitr_1.42          Formula_1.2-5      
[19] jsonlite_1.8.4      cluster_2.1.4       effectsize_0.8.3    BiocManager_1.30.20 compiler_4.2.2      httr_1.4.5         
[25] emmeans_1.8.5       backports_1.4.1     fastmap_1.1.0       lazyeval_0.2.2      cli_3.6.1           htmltools_0.5.5    
[31] tools_4.2.2         coda_0.19-4         gtable_0.3.3        glue_1.6.2          Rcpp_1.0.10         carData_3.0-5      
[37] jquerylib_0.1.4     cellranger_1.1.0    vctrs_0.6.2         nlme_3.1-160        crosstalk_1.2.0     conflicted_1.2.0   
[43] insight_0.19.1      xfun_0.38           openxlsx_4.2.5.2    timechange_0.2.0    lifecycle_1.0.3     limonaid_0.1.5     
[49] hms_1.1.3           parallel_4.2.2      yaml_2.3.7          curl_5.0.0          memoise_2.0.1       gridExtra_2.3      
[55] sass_0.4.5          rpart_4.1.19        reshape_0.8.9       stringi_1.7.12      bayestestR_0.13.1   corrplot_0.92      
[61] checkmate_2.1.0     zip_2.2.2           rlang_1.1.0         pkgconfig_2.0.3     evaluate_0.20       lattice_0.20-45    
[67] labeling_0.4.2      htmlwidgets_1.6.2   tidyselect_1.2.0    here_1.0.1          GGally_2.1.2        magrittr_2.0.3     
[73] R6_2.5.1            magick_2.7.4        generics_0.1.3      pillar_1.9.0        haven_2.5.2         foreign_0.8-83     
[79] withr_2.5.0         datawizard_0.7.1    abind_1.4-5         nnet_7.3-18         crayon_1.5.2        car_3.1-1          
[85] utf8_1.2.3          tzdb_0.3.0          rmarkdown_2.21      grid_4.2.2          readxl_1.4.2        data.table_1.14.8  
[91] digest_0.6.31       xtable_1.8-4        munsell_0.5.0       viridisLite_0.4.1   bslib_0.4.2         quadprog_1.5-8     

 

---

A work by Claudiu Papasteri

claudiu.papasteri@gmail.com